Freely available: https://earthexplorer.usgs.gov/
library(raster)
library(sp)
l.blue <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B2.TIF')
l.green <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B3.TIF')
l.red <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B4.TIF')
l.nir <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B5.TIF')
l.swir1 <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B6.TIF')
l.swir2 <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B7.TIF')
l.pc <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B8.TIF')
l.cirrus <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B9.TIF')
l.t1 <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B10.TIF')
l.t2 <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B11.TIF')
# Create raster brick
ls = brick(l.red,l.green,l.blue,l.nir,l.swir1,l.swir2,l.cirrus,l.t1,l.t2)
# CROP
e <- extent(280000, 330000, 4750000, 4800000)
ls.c = crop(ls,e)
names(ls.c) <- c('red','green','blue','NIR','SWIR1','SWIR2','cirrus','thermal1','thermal2')
writeRaster(ls.c,filename = 'Madison_Landsat_20170602.tif', format="GTiff", overwrite=TRUE)
library(raster)
## Loading required package: sp
library(sp)
r1 = brick('Madison_Landsat_20170602.tif')
r2 = brick('Madison_Landsat_20180317.tif')
names(r1) <- c('red','green','blue','NIR','SWIR1','SWIR2','cirrus','thermal1','thermal2')
names(r2) <- c('red','green','blue','NIR','SWIR1','SWIR2','cirrus','thermal1','thermal2')
r1
## class : RasterBrick
## dimensions : 1667, 1667, 2778889, 9 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 280005, 330015, 4750005, 4800015 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : /Users/hilarydugan/Documents/Zoo955/Lecture9_Landsat/Madison_Landsat_20170602.tif
## names : red, green, blue, NIR, SWIR1, SWIR2, cirrus, thermal1, thermal2
## min values : 3868, 5598, 7421, 5232, 5142, 5058, 4995, 20368, 19692
## max values : 65535, 65535, 61486, 65535, 65535, 65535, 6626, 38316, 32482
plot(r2)
plotRGB(r1, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin",
main = "Landsat True Color Composite")
plotRGB(r2, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin",
main = "Landsat True Color Composite")
# Load Lake shapefile
library(sf)
## Warning: package 'sf' was built under R version 3.4.4
## Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3
mendota <- st_read('../Lecture2_CRS/Data/LakeMendota.shp') %>%
st_transform('+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0')
## Reading layer `LakeMendota' from data source `/Users/hilarydugan/Documents/Zoo955/Lecture2_CRS/Data/LakeMendota.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 12 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -89.48369 ymin: 43.07384 xmax: -89.36737 ymax: 43.14706
## epsg (SRID): 4269
## proj4string: +proj=longlat +datum=NAD83 +no_defs
# random sample of 100 points
l.pts <- st_sample(mendota, size = 100) %>% st_sf
lake.s <- extract(r1, as(l.pts,'Spatial')) #summer values
lake.w <- extract(r2, as(l.pts,'Spatial')) #winter values
# Load watershed shapefile
watershed <- st_read('../Lecture3_Shapefiles/Data/YaharaBasins/Mendota_Basin.shp') %>%
st_transform('+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0')
## Reading layer `Mendota_Basin' from data source `/Users/hilarydugan/Documents/Zoo955/Lecture3_Shapefiles/Data/YaharaBasins/Mendota_Basin.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 552122.5 ymin: 286231.1 xmax: 584702.5 ymax: 321699.8
## epsg (SRID): NA
## proj4string: +proj=tmerc +lat_0=0 +lon_0=-90 +k=0.9996 +x_0=520000 +y_0=-4480000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
# random sample of 100 points
ws.pts <- st_sample(watershed, size = 100) %>% st_sf
ws.s <- extract(r1, as(ws.pts,'Spatial'))
ws.w <- extract(r2, as(ws.pts,'Spatial'))
# To see some of the reflectance values
lake.S <- colMeans(lake.s)
lake.W <- colMeans(lake.w)
ws.S <- colMeans(ws.s,na.rm = T)
ws.W <- colMeans(ws.w,na.rm = T)
plot(ws.S,xlab='Bands',ylab='Reflectance',type='b',xaxt='n')
axis(1,at = 1:9,labels = names(r2),las=2)
lines(lake.S,type='b',col='blue')
lines(ws.W,type='b',col='black',pch=16)
lines(lake.W,type='b',col='blue',pch=16)
legend('topleft',legend = c('Watershed - Summer','Lake - Summer','Watershed - Winter','Lake - Winter'),col = c('black','blue','black','blue'),
pch=c(21,21,16,16),bty='n')
Vegetation strongly reflects near-infrared. Plot NIR band as red
plotRGB(r1, r = 4, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "NIR Color Composite")
plotRGB(r2, r = 4, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "NIR Color Composite")
NDVI = (NIR − Red) / (NIR + Red)
ndvi1 <- (r1[[4]] - r1[[1]]) / (r1[[4]] + r1[[1]])
ndvi2 <- (r2[[4]] - r2[[1]]) / (r2[[4]] + r2[[1]])
plot(ndvi1, main = 'NDVI June')
plot(ndvi2, main = 'NDVI March')
Areas of barren rock, sand, or snow usually show very low NDVI values (for example, 0.1 or less). Sparse vegetation such as shrubs and grasslands or senescing crops may result in moderate NDVI values (approximately 0.2 to 0.5). High NDVI values (approximately 0.6 to 0.9) correspond to dense vegetation such as that found in temperate and tropical forests or crops at their peak growth stage.
As expected the NDVI is higher during the summer
Normalized Difference Water Index (NDWI), as described by McFeeters (1996), to differentiate water from non-water
NDWI <- (NIR - GREEN)) / (NIR + GREEN)
ndwi1 <- (r1[[2]] - r1[[4]]) / (r1[[4]] + r1[[2]])
ndwi2 <- (r2[[2]] - r2[[4]]) / (r2[[4]] + r2[[2]])
plot(ndwi1, main = 'NDWI June')
plot(ndwi2, main = 'NDWI March')
Normalized Difference Water Index (NDSI), to differentiate snow/ice
NDSI <- (GREEN - SWIR)) / (SWIR + GREEN)
ndsi1 <- (r1[[2]] - r1[[5]]) / (r1[[5]] + r1[[2]])
ndsi2 <- (r2[[2]] - r2[[5]]) / (r2[[5]] + r2[[2]])
plot(ndsi1, main = 'NDSI June')
plot(ndsi2, main = 'NDSI March')
Get an estimate of spatial extent of different features
NDVI: Pixels having NDVI values greater than 0.4 are definitely vegetation. Following operation masks all non-vegetation pixels.
par(mfrow=c(1,2))
ndvi1[ndvi1 < 0.4] = NA
plot(ndvi1, main = 'NDVI June')
ndvi2[ndvi2 < 0.4] = NA
plot(ndvi2, main = 'NDVI March')
NDWI: Pixels having NDVI values greater than 0.2 may be water.
par(mfrow=c(1,2))
ndwi1[ndwi1 < 0] = NA
plot(ndwi1, main = 'NDWI June')
ndwi2[ndwi2 < 0] = NA
plot(ndwi2, main = 'NDWI March')